subroutine exe_external_event(lunit,cunit,xid,pid)
  use param
  use mat_param,only:ihosts
  use global_param,only:extern_type
  use global_param,only:loop_dpa
  use global_param,only:loop_sdpa
  use box_param,only:bwidth
  use table_param,only:del_event_table
  use table_param,only:n_del_event_table
  use table_param,only:event_table
  use table_param,only:n_event_table
  use table_param,only:event_profile
  implicit none

  integer,intent(in) :: lunit,cunit,xid,pid

  integer :: i,j,idx,idc,eid,cid,sid,lid,nc,dir
  real(dbl) :: rn,rn2,re,ce,dise
  real(dbl) :: disn,dpa,casn,cas,dist
!  real(dbl) :: ox,oy,oz
  real(dbl) :: off_pos(3)
!  real(dbl) :: px,py,pz,px0,py0,pz0
  real(dbl) :: pos(3),pos0(3),vec(3),vec0(3)
!  real(dbl) :: dx,dy,dz,dx0,dy0,dz0
  real(dbl) :: d1,d2,d3,d4,d5,d6
  logical :: removed

  integer :: pos1,pos2
  character(len=len_linebuf) :: str

  if (extern_type(xid)%event_type.le.2) then
     ! <---- neutron and ion ----> !

     ! obtain the displacement energy
     dise=extern_type(xid)%dise

     ! Determine the dpa
     call random_number(rn)
     rn=rn*extern_type(xid)%dlength
     do i=1,extern_type(xid)%n_dpa
        if (rn.le.extern_type(xid)%dprob(i)) then
           re=extern_type(xid)%dgrid(i)
           idx=i
           exit
        end if
     end do

     ! Measure the recoil energy and count the event
     extern_type(xid)%dpa_e=re
     extern_type(xid)%count_dpa(idx)=&
          extern_type(xid)%count_dpa(idx)+1

     ! Calculate the dpa
     disn=0
     dpa=0
     if (re.gt.dise) then
        disn=0.8*re/(2.0*dise)
        dpa=disn*ihosts
     end if
     extern_type(xid)%partial_dpa=dpa
     loop_dpa=dpa

     ! Determine the cascade sample
     call random_number(rn)
     rn=rn*extern_type(xid)%clength
     do i=1,extern_type(xid)%n_cas
        if (rn.le.extern_type(xid)%cprob(i)) then
           ce=extern_type(xid)%cgrid(i)
           idc=i
           exit
        end if
     end do

     ! Measure the cascade energy and count the event
     extern_type(xid)%sub_e=ce
     extern_type(xid)%count_cascade(idc)=&
          extern_type(xid)%count_cascade(idc)+1

     ! Determine the cascade energy
     casn=0
     cas=0
     if (ce.ge.dise) then
        casn=0.8*ce/(2.0*dise)
        cas=casn*ihosts
     end if
     extern_type(xid)%partial_sdpa=cas
     loop_sdpa=cas

     ! Do nothing 
     ! if the cascade energy is less than displacement energy
     if (ce.lt.dise) return

     ! Do nothing
     ! if the cascade samples does not exist.
     if (extern_type(xid)%n_objset.eq.0) return

     ! Do nothing
     ! if the cascade samples does not exist in this group.
     nc=extern_type(xid)%cascade_list(idc)%nlist
     if (nc.eq.0) return

     ! Choose the cascade sample among the list in this group
     call random_number(rn)
     cid=ceiling(rn*nc)
     sid=extern_type(xid)%cascade_list(idc)%clist(cid)

#ifdef DEBUG
     pos1=0
     pos2=pos1+19
     write(str(pos1+1:pos2),&
          '("cascade set no.=",i3)') cid
     pos1=pos2
     call log_msg1(cunit,"")
     call log_msg1_indent(cunit,1,str(1:pos2))
#endif

     select case(extern_type(xid)%rand)
     case(0)
        ! Put the cascade sample without manipulation.
        do i=1,extern_type(xid)%objset(sid)%n_obj
           lid=extern_type(xid)%objset(sid)%obj(i)%obj_type
           do j=1,3
              pos0(j)=extern_type(xid)%objset(sid)%obj(i)%pos(j)
           end do
           dir=extern_type(xid)%objset(sid)%obj(i)%dir
           do j=1,3
              vec(j)=extern_type(xid)%objset(sid)%obj(i)%vec(j)
           end do
           call boundary_check(pos0(1),pos0(2),pos0(3),removed)
           if (removed) cycle
           call create_object(lid,pos0(1),pos0(2),pos0(3),&
                dir,vec(1),vec(2),vec(3))
        end do
     case(1)
        ! Put the cascade sample with random offsets
        ! without rotation
        call random_number(rn)
        off_pos(1)=rn*bwidth(1)
        call random_number(rn)
        off_pos(2)=rn*bwidth(2)
        call random_number(rn)
        off_pos(3)=rn*bwidth(3)
        do i=1,extern_type(xid)%objset(sid)%n_obj
           lid=extern_type(xid)%objset(sid)%obj(i)%obj_type
           do j=1,3
              pos(j)=off_pos(j)+&
                   extern_type(xid)%objset(sid)%obj(i)%pos(j)
           end do
           dir=extern_type(xid)%objset(sid)%obj(i)%dir
           do j=1,3
              vec(j)=extern_type(xid)%objset(sid)%obj(i)%vec(j)
           end do
           call boundary_check(pos(1),pos(2),pos(3),removed)
           if (removed) cycle
           call create_object(lid,pos(1),pos(2),pos(3),&
                dir,vec(1),vec(2),vec(3))
        end do
     case(2)
        ! Put the cascade sample with random rotation
        ! determine the offset
        call random_number(rn)
        off_pos(1)=rn*bwidth(1)
        call random_number(rn)
        off_pos(2)=rn*bwidth(2)
        call random_number(rn)
        off_pos(3)=rn*bwidth(3)

        ! alpha
        call random_number(rn)
        d1=cos(2.0*rn*pi)
        d2=sin(2.0*rn*pi)
        ! beta
        call random_number(rn)
        d3=cos(2.0*rn*pi)
        d4=sin(2.0*rn*pi)
        ! gamma
        call random_number(rn)
        d5=cos(2.0*rn*pi)
        d6=sin(2.0*rn*pi)

        do i=1,extern_type(xid)%objset(sid)%n_obj

           lid=extern_type(xid)%objset(sid)%obj(i)%obj_type

           ! original position
           do j=1,3
              pos0(j)=extern_type(xid)%objset(sid)%obj(i)%pos(j)
           end do

           ! rotate in x-axis
           pos(1)=pos0(1)
           pos(2)=pos0(2)*d1-pos0(3)*d2
           pos(3)=pos0(2)*d2+pos0(3)*d1

           ! rotate in y-axis
           pos0(1)=pos(1)*d3+pos(3)*d4
           pos0(2)=pos(2)
           pos0(3)=-pos(1)*d4+pos(3)*d3

           ! rotate in z-axis
           pos(1)=pos0(1)*d5-pos0(2)*d6
           pos(2)=pos0(1)*d6+pos0(2)*d5
           pos(3)=pos0(3)

           ! apply the offsets
           do j=1,3
              pos(j)=pos(j)+off_pos(j)
           end do

           call boundary_check(pos(1),pos(2),pos(3),removed)

           if (removed) cycle
           
           dir=1
           do j=1,3
              vec(j)=0
           end do

           ! treat the direction
           if (extern_type(xid)%objset(sid)%obj(i)%dir.eq.0) then
              ! The direction is given.
              dir=0

              ! original position
              do j=1,3
                 vec0(j)=extern_type(xid)%objset(sid)%obj(i)%vec(j)
              end do

              ! rotate in x-axis
              vec(1)=vec0(1)
              vec(2)=vec0(2)*d1-vec0(3)*d2
              vec(3)=vec0(2)*d2+vec0(3)*d1

              ! rotate in y-axis
              vec0(1)=vec(1)*d3+vec(3)*d4
              vec0(2)=vec(2)
              vec0(3)=-vec(1)*d4+vec(3)*d3

              ! rotate in z-axis
              vec(1)=vec0(1)*d5-vec0(2)*d6
              vec(2)=vec0(1)*d6+vec0(2)*d5
              vec(3)=vec0(3)

              do j=1,3
                 if (vec(j).gt.0) then
                    vec(j)=1
                 else
                    vec(j)=-1
                 end if
              end do

           end if
           ! create the object
           call create_object(lid,pos(1),pos(2),pos(3),&
                dir,vec(1),vec(2),vec(3))
        end do
     case default
     end select

     ! Count the occurence of this event
     nc=extern_type(xid)%n_occur
     nc=nc+1
     extern_type(xid)%n_occur=nc
     if (nc.ge.extern_type(xid)%max_n_occur) then
        ! remove the event from the table
        n_del_event_table=n_del_event_table+1
        eid=event_profile(pid)%event_list(1)
        del_event_table(n_del_event_table)=eid
        event_table(eid)%removed=.true.
        n_event_table=n_event_table-1
        event_profile(pid)%n_event=0
     end if

  end if

  if (extern_type(xid)%event_type.eq.3) then
     ! <---- electron ----> !

     ! obtain the displacement energy
     dise=extern_type(xid)%dise

     ! determine the dpa
     call random_number(rn)
     rn=rn*extern_type(xid)%dlength
     do i=1,extern_type(xid)%n_dpa
        if (rn.le.extern_type(xid)%dprob(i)) then
           re=extern_type(xid)%dgrid(i)
           idx=i
           exit
        end if
     end do

     ! Measure the recoil energy and count the event
     extern_type(xid)%dpa_e=re
     extern_type(xid)%count_dpa(idx)=&
          extern_type(xid)%count_dpa(idx)+1
     
     ! Calculate the dpa
     disn=0
     dpa=0
     if (re.gt.dise) then
        disn=0.8*re/(2.0*dise)
        dpa=disn*ihosts
     end if
     extern_type(xid)%partial_dpa=dpa
     extern_type(xid)%partial_sdpa=dpa
     loop_dpa=dpa
     loop_sdpa=dpa

     if (re.lt.dise) return

     ! determine the default position
     call random_number(rn)
     off_pos(1)=rn*bwidth(1)
     call random_number(rn)
     off_pos(2)=rn*bwidth(2)
     call random_number(rn)
     off_pos(3)=rn*bwidth(3)

     call boundary_check(off_pos(1),off_pos(2),off_pos(3),removed)

     if (.not.removed) then
        lid=extern_type(xid)%source_type(1)
        dir=1
        do i=1,3
           vec(i)=0
        end do
        call create_object(lid,off_pos(1),off_pos(2),off_pos(3),&
             dir,vec(1),vec(2),vec(3))
     end if

     do i=2,extern_type(xid)%n_type

        if (extern_type(xid)%dist(1).gt.&
             extern_type(xid)%dist(2)) then
           call random_number(rn)
           rn=rn*(extern_type(xid)%dist(1)-&
                extern_type(xid)%dist(2))
           dist=rn+extern_type(xid)%dist(2)
        else if (extern_type(xid)%dist(1).lt.&
             extern_type(xid)%dist(2)) then
           call random_number(rn)
           rn=rn*(extern_type(xid)%dist(2)-&
                extern_type(xid)%dist(1))
           dist=rn+extern_type(xid)%dist(1)
        else
           dist=extern_type(xid)%dist(1)
        end if

        call random_number(rn)
        rn=rn*2.0*pi
        call random_number(rn2)
        rn2=rn2*2.0*pi
        d1=cos(rn2)

        pos(1)=off_pos(1)+dist*cos(rn)*d1
        pos(2)=off_pos(2)+dist*sin(rn)*d1
        pos(3)=off_pos(3)+dist*sin(rn2)

        call boundary_check(pos(1),pos(2),pos(3),removed)

        if (removed) cycle

        lid=extern_type(xid)%source_type(i)

        dir=1
        do j=1,3
           vec(j)=0
        end do

        call create_object(lid,pos(1),pos(2),pos(3),&
             dir,vec(1),vec(2),vec(3))

     end do

     ! Count the occurence of this event
     nc=extern_type(xid)%n_occur
     nc=nc+1
     extern_type(xid)%n_occur=nc
     if (nc.ge.extern_type(xid)%max_n_occur) then
        ! remove the event from the table
        n_del_event_table=n_del_event_table+1
        eid=event_profile(pid)%event_list(1)
        del_event_table(n_del_event_table)=eid
        event_table(eid)%removed=.true.
        n_event_table=n_event_table-1
        event_profile(pid)%n_event=0
     end if

  end if

end subroutine exe_external_event
